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Abstract We address the problem of stability of mo¬ 
tor actions implemented by the central nervous system 
based on simple algorithms potentially reflecting physi¬ 
cal (including physiological) processes within the body. 
A number of conceptually simple algorithms that solve 
motor tasks with a high probability of success may be 
based on feedback schemes that ensure stability of sub¬ 
spaces of neural variables associated with accomplish¬ 
ing those tasks. The task is formulated in terms of linear 
constrains imposed either on the human body mechan¬ 
ical variables or on neural variables; we discuss three 
reference frames relevant to these processes. We discuss 
underlying basic principles of such algorithms, their ar¬ 
chitecture, and efficiency, and compare the outcomes of 
implementation of such algorithms with the results of 
experiments performed on the human hand. 


Contents 


1 Introduction: The Nervous System, Motor Control, 

and Stability Search Problem. ^ 

2 The Stability Search Algorithms, Dynamics Con¬ 
trol, and the Effect of Noise. |3] 

3 Experiments with Human Hand. 1131 


4 Discussion: The Stability Search Algorithm Hypothe¬ 
ses within the Context of Motor Control Hypotheses 1181 


1 Introduction: The Nervous System, Motor 
Control, and Stability Search Problem 

Many functions of the central nervous system (CNS) 
can be described as combining numerous elements (we 
will refer to their outputs as elemental variables) into 
relatively low-dimensional sets related to such functions 
as cognition, perception, and action. The existence of 
such low-dimensional sets ensures stability of percepts, 
thoughts, and actions despite the variable contributions 
from the elements (sensory receptors, neurons, motor 
units, etc.) and changes in the environment. Here, we 
try to offer a mathematical description of processes that 
could bring about such stability based on variability. 
The aim of this paper is not to mimic a control algo¬ 
rithm for some particular task, which the CNS solves, 
but rather to find a conceptual framework, though vague 
at the moment, which would allow one to offer mathe¬ 
matical principles that can shape neural activity associ¬ 
ated with a variety of tasks. We use, as an example, the 
production of voluntary movements by redundant sets 
of elements. Our approach is not based on an a priori 
cost function, but on an intuitive, simple algorithmic 
principle. Of course, starting from a very complex al¬ 
gorithm may have an a priori advantage in achieving 
high probability of success as compared to simpler al¬ 
gorithms. We start building a model from a linear dy¬ 
namic system and then add operational complexity as 
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needed to achieve a reasonable success probability. At 
any level of description, the system for movement pro¬ 
duction is apparently redundant [Jj. This means that 
the number of its elemental state variables is larger 
than the number of constraints associated with typ¬ 
ical tasks. Depending on the selected level of analy¬ 
sis, elements and elemental variables are defined differ¬ 
ently (reviewed in [IB]). For example, kinematic analy¬ 
sis of multi-joint action frequently considers individual 
joint rotations as elemental variables. Kinetic analysis 
of the force and moment production by parallel chains 
(such as the digits of the human hand) may consider 
forces/moments produced by the digits as elemental 
variables. Analysis of muscle activation may be based 
on firing patterns of individual motor units as elemental 
variables. And so on. The apparent redundancy typical 
of all these examples has been re-cast as abundance m 
HU to emphasize that the extra elemental variables are 
not the sources of computational problems for the CNS 
but an essential component of the design that allows 
combining stability of actions with flexibility of perfor¬ 
mance (changing actions, responding to perturbations, 
performing several actions in parallel, etc.). 

Consider the following example as an illustrations 
of our goal and approach. Imagine a walking person 
who suddenly steps on a slippery surface.The slip is 
typically followed by a very complex pattern of move¬ 
ments of all body parts resulting in restoring balance 
in a large percentage of cases. Each time a slip oc¬ 
curs the movement pattern looks unique. We assume 
here that such highly variable patterns emerge as a re¬ 
sult of a single, relatively simple algorithm applied to 
cases with varying initial conditions at the slip. What 
could be the search algorithm for new stability, presum¬ 
ably formulated as a set of simple rules realized by the 
physical/physiological system without explicit compu¬ 
tational steps ? What is the success rate of not falling 
yielded by such an algorithm ? We expect a balance be¬ 
tween the operational simplicity of the algorithm and 
the success probability. What could be the relation of 
such an algorithm to established notions in the field of 
motor control such as hierarchical control and uncon¬ 
trolled manifold hypothesis [27] ? 

In the most simple mathematical setting, the equi¬ 
librium control problem (task) can be seen as a linear 
requirement Yj = a j,i x i imposed to N elemental 

variables Xi (with i running from 1 to N) describing 
states of a redundant set of effectors, which depends on 
desired position of the object given by a set of variables 
Yj (with j running from 1 to M). The task is to make 
the subspace determined by this set of linear conditions 
mechanically stable, that is to find a proper feedback 


matrix Pi.j, which via the equation 

M 

x i = J2Pi,j(yj- Y j) (!) 

t=i 

relates the time derivatives Xi of the elemental variables 
and the deviations Syj = yj — Yj of instantaneous co¬ 
ordinates yt of the controlled object from their desired 
positions Yj. One has to construct a simple algorithm 
of yielding such a feedback matrix /3ij for an arbitrarily 
predetermined task matrix ctj.i- 

If now one eliminates the physical aspect of the 
problem and treats the variables Xi and yj as certain 
brain variables responsible for the control of a physical 
object, then the model begins to describe a purely men¬ 
tal process that operates with some neural variables Xi 
responsible for the effector activity and the neural vari¬ 
ables yj responsible for the perception of the controlled 
object state, whereas the matrix ay,; connects these two 
groups of variables. Here, we suggest a sequence of gen¬ 
eral principles that may be seen as extensions of the idea 
that biological systems are reasonably sloppy DEM!- 
Under this expression, we mean that approximate solu¬ 
tions may be generated for typical tasks as long as they 
lead to success with acceptable probability. This no¬ 
tion is closely related to the notion of satisficing coined 
by Herbert Simon [29]. This means that they do not 
solve problems exactly but rather use simplified rules 
that produce solutions that are good enough (e.g., suc¬ 
cessful most of the time). From the mathematical point 
of view, even for a simple linear model, construction 
of a feedback algorithm is a nonlinear process, which 
is meant to vary the matrix elements of 8i,j depend¬ 
ing on the size and the evolution history of variables 
Xi . First, we suggest a requirement, that the algorithm 
should rely, as much as possible, on local actions, that 
is, on actions that change the way a neural variable Xi 
participates in the feedback, such that the decision to 
undertake an action is based on the actual and previ¬ 
ous values Xi(t < t actua i) of this very variable, whereas 
other variable do not affect the decision. However, such 
a purely local algorithm may be incapable of providing 
good stability properties, and thus one needs to involve 
non-local actions, when the participation of one vari¬ 
able depends on the state and the evolution history of 
other variables. In such a case, our second principle is 
that the number of nonlocal interventions should be 
kept at a minimum. Third, we would like to formulate 
the decision criteria in simple, intuitive terms. In this 
study, we start with a rule ’’Act on the most nimble” 
(the AMN-rule), when changes in the local parameters 
occur for the most unstable variable first. 
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Two more points should be discussed in the con¬ 
text of the requirements imposed on the algorithm, (i) 
It is desirable that the algorithm ensures robustness 
of the feedback and yields such that can be used 
across the task a h i within certain limits. Formally this 
implies that the spectrum {cu m } of eigenvalues of the 
matrix ojjj> = 'Y^i =1 remains stable (with nega¬ 

tive real parts) even when the task matrix experiences 
a perturbation Sctj.i. (ii) It is hard to imagine a direct 
experimental access to the assumed neural variables. 
Since the classical works by the group of Georgopoulos 
HI, studies of cortical neuronal population vectors 
have provided plenty of evidence that such vectors cor¬ 
relate with kinematic movement characteristics 121331 - 
We would like to emphasize, however, that such corre¬ 
lations are strongly dependent on the external condi¬ 
tions of the experiments and do not prove the physical 
meaning of the population vectors. In contrast, we are 
implying neural variables directly related to parameters 
related to setting a task in a multi-dimensional space 
of elemental variables and participating in stabilization 
of its performance. Fortunately, analysis of the mathe¬ 
matical structure of the problem can rely on the noise 
statistics. The reason for this is the fact that, for a 
rather general model of noise, one can find a Gaussian 
distribution of the measurements in the M -dimensional 
space of the observed variables yj, while the eigenvec¬ 
tors of the covariance matrix of this distribution coin¬ 
cides with the eigenvectors of uijj >, and a certain rela¬ 
tion among the eigenvalues can be deduced. 

The ” act on the most nimble” rule makes our model 
essentially non-linear. However, in-between successive 
actions on the most nimble variables, the dynamical 
system can be treated as linear and open, which im¬ 
plies its being subjected to noise, defined as uncontrol¬ 
lable small and random deviations of the coordinates 
from their averages over short time intervals. Though 
biological systems typically do not manifest Gaussian 
noise statistics, the basic characteristics, such as noise 
covariance matrix, still can reveal important instanta¬ 
neous characteristics of the matrix governing the linear 
part of dynamics, such as eigenvectors and, to a certain 
extent, eigenvalues, thus helping to compare calculated 
and observed results. 

2 The Stability Search Algorithms, Dynamics 
Control, and the Effect of Noise 

In this Section we discuss the main ideas underlying 
construction of the stability search algorithms. We start 
by discussing the mathematical complexity of the prob¬ 
lem, and specific requirements for solving this problem 
by an analog computer - a model for the CNS at this 


stage. In this context we talk about local algorithms 
contrasting them to the non-local ones and introduce 
a new basis of neural variables where such local algo¬ 
rithms can be implemented. 

We continue by exploring the ingredients required 
for the stabilization algorithm construction. By propos¬ 
ing a simple model of the control we specify the prob¬ 
lem in terms of a task and a feedback matrices, showing 
that a random choice of the feedback has extremely low 
chances of yielding the required stabilization. We next 
turn to the simplest local equilibrium search algorithm 
- ” act on the most nimble”, which can considerably in¬ 
crease the chances of stabilization gaining. This idea 
turns to be even more efficient when the feedback ma¬ 
trix is not random, but is tailored in a way which we call 
’’generations” structure, - where each one-dimensional 
subspace is coupled to a multidimensional ” generation” 
of variables, that all have typical velocities (’’nimble¬ 
ness”) of the same order of magnitude. Each ’’genera¬ 
tion” has ’’nimbleness” considerably slower that for the 
previous one and much higher than for the next ” gener¬ 
ation” . This idea is further developed when we account 
for the possibility of non-local actions, where each next 
generation, when started at a relevant time scale to im¬ 
plement local operations, is also allowed to inhibits ac¬ 
tivity of previous generations. This is followed by anal¬ 
ysis of the case where each generation nimbleness is 
a tunable fit parameter, and a solution can always be 
found, for a system of any dimensionality, although this 
solution rely on the typical velocities that may differ by 
many orders of magnitude for the first and the last gen¬ 
erations. 

The next topic we discuss in this Section is the re¬ 
lations between noise and stability, which later on will 
help us to gain a deeper insight into the dynamics of 
the system by analyzing experimentally observed co- 
variances of the noise. The term “noise”is ambiguous in 
neurophysiology (reviewed in [35] . For example, during 
quiet standing, humans show spontaneous deviations of 
the body from the vertical (postural sway), which are a 
superposition of two components, rambling and trem¬ 
bling [23] . The former reflects migration of the equilib¬ 
rium point with respect to which balance is organized, 
while the latter reflects oscillations about that moving 
equilibrium point. Rambling has been viewed as a re¬ 
flection of an unintentional but meaningful search for 
limits of stability and, hence, it would be counterin¬ 
tuitive to classify it as “noise”. Nevertheless, for the 
purposes of the current study, we will use “noise” to im¬ 
ply spontaneous variance in physiological signals that 
is not related to the explicit task and is not induced by 
an identifiable external perturbation. We understand 





4 


Please give a shorter version with: \authorrunning and \titlerunning prior to \maketitle 


that this definition may attribute to noise physiologi¬ 
cally meaningful processes. 

We conclude the Section by considering a realistic 
realization of the control algorithm in hierarchical sys¬ 
tems, show that for this architecture of the feedback 
damping of the dynamic variables is indispensable and 
propose a model of self-adjusting damping, which yields 
time dependencies of the dynamic variables resembling 
ones observed experimentally. 


2.1 Control Complexity and Coordinate Frames 

Finding solution of M linear equations Yj = Y^iLi a j,i x i 
for N > M variables is a problem of a polynomial 
complexity ~ M 2 . It is sufficient to take the first M 
columns of the M x N matrix oy 7 and invert the re¬ 
sulting M x M square matrix, provided it is not degen¬ 
erate. Solutions of the problem for the first M variables 
{x\... xm} form a M-dimensional subspace Sm in the 
entire iV-dimensional space, parametrized by the re¬ 
maining N — M variables {xm+i ■ ■ ■ %n}- Construction 
of the inverse matrix a~j based on the sequential find¬ 
ing of its line vectors complemented by the orthonor¬ 
malization of these vectors with those found earlier in¬ 
deed implies the number of operations on the order of 
M 2 . 

To solve this problem on an analog computer, one 
should construct and implement a dynamic process that 
has the M-dimensional subspace Sm as a stable station¬ 
ary manifold, as suggested by Eq. 0. Apparently, the 
additional requirement of stability augments the math¬ 
ematical complexity of the relevant algorithm, but once 
this requirement is fulfilled, the result offers an advanta¬ 
geous generalization - it can be implemented for solving 
systems of nonlinear equations that is a problem of a 
much higher complexity. Here, we are going to propose 
several strategies of constructing heuristic algorithms 
that result in the dynamic stabilization of subspaces 
defined by systems of M equations for N > M vari¬ 
ables. 

Two points have to be elucidated. The first is that, 
according to the general vision, an analog computer is 
a collection of rather simple, although linearly interact¬ 
ing, dynamic systems - elements, while the algorithms 
are viewed in this context as sequences of instrumental 
prescriptions of how to modify these interactions. Each 
of the prescribed conditions is imposed on the value 
of a dynamic variable describing the state of one ele¬ 
ment and the prescribed action is applied to the same 
or to another element. Prescribed action on one chosen 
element that depends only on the state of either that 
element or on another chosen element implies so-called 
locality or bi-locality, respectively, of the elementary 


step of the algorithms. The second point is that non¬ 
linearity of the analog computer generated by such a 
local and/or bi-local algorithm may not result in univer¬ 
sal stability of the dynamics. In other words, for some 
initial conditions the system may reach a stationary 
solution, while for other initial conditions it becomes 
unstable. In such a situation, one can speak about the 
success rate R(N , M) of the algorithm, and this statis¬ 
tical quantity serves as its quality characteristic. 

We now move to implementing the two principles 
mentioned in the Introduction: (1) Reliance primarily 
on local actions while nonlocal actions are minimally in¬ 
volved; and (2) Using the intuitive AMN-rule. First, we 
introduce three types of specific bases related to certain 
linear combinations of the body variables. First, there 
is a measurement basis formed by experimentally acces¬ 
sible variables, for example positions, forces, or muscle 
activations. Second, there is another, task-dependent, 
basis formed by so-called modes , that are linear com¬ 
binations of the former variables independently fluc¬ 
tuating under the action of noise and coinciding with 
the eigenvectors of the aforesaid matrix Ujj> (cf. [filll4|. 
The modes are expected to be stable for broad ranges 
of tasks and not to change without specialized training. 
The third basis is the reference system relying on the 
variables Xi responsible for the control over the body 
state. To our knowledge, this coordinate system has 
not been defined previously. Each of the coordinates 
of the third system represents combinations of modes 
that are task specific and relatively quickly adjustable 
to changes in the external conditions of task execu¬ 
tion, for example to changes in stability requirements 
(e.g. 0). This basis comprises the variables that ex¬ 
perience local control, that is, changes in the feedback 
loop gain with respect to each of these variables Xi pro¬ 
duced once a certain functional Zi{x.i(t)}oi these very 
variable reaches a critical value. Such gain changes may 
occur gradually with changes in the value of the func¬ 
tional. To summarize, the task is considered within a 
preexisting, task-dependent basis. 


2.2 The Control Structure 


To solve the problem of a subspace stabilization one 
first takes an arbitrarily M x N matrix a that imposes 
conditions of Eq.|l]) with Yj = 0 thus defining the re¬ 
quired subspace Sm = {Vi = 0}. Thus the task is to find 
a, N x M linear feedback matrix /3 that, by determining 
the derivatives 


M 


d 
dt 

3 =1 


PijVji 


( 2 ) 
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prescribes the modifications of {xi} leading towards the 
subspace Sm- Dynamics of the system is therefore given 
by the equation 


d 

dt 


Xi = 


N M 

EE 

3 =1 k =1 


(3) 


and a stable M -dimensional subspace Sm exists in the 
TV-dimensional space if all the non-zero eigenvalues of 
the degenerate square TV x TV matrix Q = /3a of the 
rank M have negative real parts. The probability to 
satisfy this requirement for randomly chosen either a 
or both a and /3 is rather low, and in order to stabi¬ 
lize Sm with a high enough success rate R(N , M) for 
a chosen a, a randomly or systematically chosen lin¬ 
ear feedback matrix initial may need to be modified 
Anitiai /3 m odified once it does not ensure negativity of 
the /3a eigenvalues real parts and hence does not yield 
the required stabilization. A sequence of the modifi¬ 
cations /3 mo dified(n) t /3 mo dified(n+i) can be continued 
according to a certain algorithm until stabilization is 
achieved. 

Each elementary act of the modification algorithm 
rely on three main ingredients: (i) the variable x j initiat¬ 
ing the action, (ii) type of the action itself /3 mo difl e d(n) = 
G/3 mo dified(n+i) implementing the modification A„odified(n 
/3 m odified(n+i), which is specified by a non-linear func¬ 
tional matrix G({x/(f)} , {^(f)}), and (iii) the variable 
Xf subject to the action. Once the source and the tar¬ 
get variables coincide, we encounter a so-called local 
algorithm, which implies that for f ^ i we deal with 
a non-local elementary act of the algorithm. During 
last decades, the difference between local and non-local 
actions are has been widely discussed in the context 
of Quantum Informatics, where physical implementa¬ 
tion of non-local transformations, so-called quantum 
gates, is a much more challenging experimental task 
as compared to local transformations. In a situation 
where each variable Xk is associated with some location 
in space, the situation typical of Quantum Informat¬ 
ics becomes generic: Local algorithms are much easier 
to implement in a physical realization of the controlled 
systems than non-local ones since they do not require 
transportation of the nonlinear action from one spatial 
location to another. 


will thus be employed as entries into the feedback loop 
given by a rectangular TV x M matrix ffj, such that 
Eq. © holds, while Eq. © describes the system dynam¬ 
ics. 

For a generic rectangular randomly chosen matrix 
a.ij and the randomly chosen feedback matrix /3j ; fc the 
success rate i?(TV, M) found numerically scales approxi¬ 
mately as 2 -M , which is consistent with the assumption 
that all non-zero eigenvalues of the matrix 17 may have 
the real parts both negative and positive with equal 
probability and, roughly speaking, different eigenvalues 
are statistically independent one from one another. 


2.2.2 Local algorithm for random feedback 


Our aim now is to construct a control algorithm, which 
is capable of improving the success rate. We call this 
algorithm the AMN-rule later. We begin with the sim¬ 
plest case of a local algorithm by assuming that the 
feedback sign for a given x^ changes once a positive 
quantity 

a(i) = i (t 1 ) dt (4) 

exceeds a threshold value Z. This represent an example 
) “df ’’act on the most nimble” control, where, as a first 
guess, we consider the integral of the rates squared, 
which is an analog of the strictly positive dissipated 
heat functional typical of electrical circuits. Further, 
we will consider other strictly positive functionals. For¬ 
mally, the local control implies that the set of equations 
([3]) is modified by the presence of a sign function and 
reads 

d N M 

jXi = G u EE Pi,k^k,jXj (5) 

0 = 1 k —1 

Gu = sign(Z - Zi(t)). (6) 


This strategy yields 100% success rate for the case 
of M = 1, that is the case where the matrix /3a has rank 
1. In this case the dynamic process occurs in a subspace 
of the dimension 1 such that the dynamic equation 



N 

oii t iGupi t iyi 

i=l 


(7) 


2.2.1 No algorithm 

We first consider the simplest realization of an analog 
computer working without a controlling algorithm. For 
the sake of presentation simplicity, assume that the de¬ 
sired subspace Sm corresponds to zero values of all the 
variables //ie{i,M}- Non-zero values of these variables 


for the variable y\ = XEi a i,i x i determines time de¬ 
pendencies of all the other variables Xi = In¬ 
deed, if the quantity ui = is positive, jq 

exponentially increases in time. This means that the 
variable Xi corresponding to the largest produces 
the highest ’’Joule’s heat” z-ft) and at some moment 
of time, when this quantity exceeds the threshold value 
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M 

1 

2 3 

4 

4 

4 

5 

N 

10 

20 30 

15 

40 

60 

30 

^random 

0.50 0.24 0.14 

- 

0.06 

- 

- 

-^local 

1 

0.77 0.46 

0.17 

0.16 0.154 

0.06 

-^tailored 

1 

0.85 0.62 0.32/0.28 0.36 

0.31 

0.25/0.21 

-^generations 

1 

0.9 0.68 0.42/0.46 0.66 

- 

0.45/0.31 


Table 1 The success rate of various algorithms stabilizing M-dimensional subspace of IV-dimensional space. 


Z, the sign of sign (Z — Zi(t))Bi,i changes. For a positive 
this results in a decrease of w and slowing-down 
of the instability. For a random matrix a, however, the 
matrix element aq^ can equally be negative or posi¬ 
tive; in the latter case, the change of the sign results in 
the instability speeding-up. In this case, the exponential 
growth of the variables Xi continues, and after awhile, 
the next biggest \Bi,i\ leads to a change of the corre¬ 
sponding feedback sign. Changes of the signs continue 
till u) becomes negative, and the dynamic process be¬ 
comes stable. This will definitely occur for the matrix /? 
of the rank 1, but for higher ranks this stabilization may 
not happen. Moreover, numerical search shows that the 
success rate i?i OC ai of such stabilization drops with the 
increasing rank M. The results of the numerical sim¬ 
ulation for this case are shown in Table [T] the fourth 
row (i?iocai) ■ These results have to be compared with 
ifrandom for the not controlled random feedback. While 
the success rate is higher for the controlled system for 
all M and N, it drops rather quickly with the rank M, 
while being relatively less sensitive to N (compare the 
R values for M = 4 and N = 15, 40, and 60). 

2.2.3 Tailored linear feedback and the local algorithm 

Remaining within the general case of Eq. ([5]) , we replace 
the random feedback matrix Bi.j by another one, con¬ 
structed to augment the success rate for higher M. We 
take it in the form 




0 

0 . 

' ° ^ 



0 

qrf 

0 . 

0 


B = 

0 

0 

q 2 lf . 

0 

(8) 


V 0 

0 

0 . 

. ) 



where it is a column vector of the size of the ratio 
N/M , which we assume to be an integer. Here q is a 
small parameter, which ensures that at each time scale 
~ q~ n , one has to deal with a subspace of rank just 
1 with the same algorithm as for the case of M = 1, 
while the subspaces corresponding to the larger matrix 
elements ~ q~ k<n are assumed to have already been 
stabilized earlier. 


In simple words the main idea is: there are genera¬ 
tions of elemental variables organized by their expected 
nimbleness, which scales as q 9 , that is, exponentially 
with the number g of the generation. The controller uses 
the AMN-rule within an appropriate channel (which is 
a set of variables aq grouped by their common input 
from a component of the vector l/) and then deals with 
the next most nimble channel without returning to the 
first one, and so on. This approach improves the success 
rate for higher M ; however, even it cannot guarantee 
100% convergence. The results of the numerical work 
for the success rate -Rtaiiored in the case M = 3, N = 9, 
q = 0.07, 



are shown on the fifth line of Table [lj while the results 
for q = 0.1 calculated for two of these cases are shown 
after the slash. 


2.2.4 Tailored linear feedback and the 
generation-specific non-local algorithm 

The idea of bi-local control allows excluding undesired 
changes in the feedback sign determined at an earlier 
stage of control that may be induced at a later stage. 
In the feedback matrix Eq.(|8|, we identify parts that 
belong to different generations, corresponding to dif¬ 
ferent orders of the parameter q. The first generation 
corresponds to q°, and the last, the most recent, gen¬ 
eration to q AI ~ l . Evidently, each generation accounts 
for the feedback at the corresponding time scale. The 
idea of the control algorithm is that changing the sign 
of a variable belonging to a generation blocks changes 
of the signs of the variables belonging to all former gen¬ 
erations. 

Formally, the bi-local control implies that the set 
of equations ( |4|5[ ) is modified by the presence of step 
factors 0(Z — zf) at the derivatives for the variables 
Zi corresponding to ’’generations” that happened after 
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than that of i. and reads 


2.3 Feedback Strength Exploring Algorithm 


d_ 

dt 


Zi(t) 


dt Xi 


G, 


II 9{Z-Zr{t)) 

rG inferior (i) 

N M 

Gu EE (3i,k&k,j%j 

j =1 k =1 
sign(Z - 



( 10 ) 


where 0{x) is the step function. The first equation of 
this system shows that the functionals Zi(t) governing 
the feedback signs are no longer local, since dynamic 
equations ruling these quantities depend not only on the 
corresponding local squared velocities but also on the 
values of the functionals for other variables. This con¬ 
struction further improves the success rate. The results 
of the numerical work for the success rate -R ge nerations 
in this case are shown in the last line of the Table Q] for 
q = 0.07-0.1. 


y(t) 



Fig. 1 An example of stabilizing a four-dimensional sub¬ 
space in a 40-dimensional space. Time dependencies of the 
parameters y 1 ... 4 , that determine the subspace for yi = 0 are 
shown. The first variable y\ (the blue solid line) corresponds 
to the first generation of the variables xi ... no having the 
strongest coupling. The second variable, j/ 2 , (the red dashed 
line) is coupled to the variables m • • - no with weaker cou¬ 
pling constants, scaled by the factor q = 0.4 with respect to 
the first variable. The third and the fourth generations (the 
dash-dot and the dotted lines, green and black, respectively) 
have couplings scaled by the factors q 2 and <j 3 , respectively. 
In the figure, one can see seven sequential discontinuities of 
the derivatives of the dependencies related to the changes of 
signs of the most nimble variables Xi , that first occurs in the 
first generation, then in the second, etc. The changes in each 
of the generations manifest themselves in all dependencies via 
the matrix a 


In Figf^we show an example of the time dependent 
deviations yi... 4 (t) in the case of a successful control for 
N = 4, M = 40, q = 0.4, where the abrupt changes of 
the dependencies reflect changes of the feedback signs. 


The bi-local algorithm can be modified to gain the 100% 
success rate if instead of the fixed power factors q n en¬ 
tering Eq.([8]) one allows sequential choosing of these 
feedback strength parameters for each next generation. 
I 11 a sense, this algorithm implies learning, that is, given 
a task a , it modifies the feedback matrix /3 once the 
sign-changing algorithms does not result in the required 
subspace stabilization. More specifically, given q < 1, 
similar to Eq.([8| where 


P = 


/it 0 0 

0 q Pl lt 0 
0 0 q P2 it 


° \ 
0 


V 0 0 


q PM-iit J 


( 11 ) 


the power parameters pi < P 2 < .. - Pm- 1 that have to 
be chosen such that the subspace {?/» = 0} is stable. 

One begins with all pi = 00 , that is, with a one¬ 
dimensional subspace, which can always be made sta¬ 
ble by proper choice of signs. Next, one sets pi to zero, 
and starts to implement the sign changing algorithm 
in the subspace of the second diagonal cell of Eq. (13). 
If this algorithm does not lead to a stable subspace of 
the dimension 2, one increases p\ by unity, and imple¬ 
ments the sign-changing algorithm for the second cell 
once again. Repeating this procedure leads to finding 
Pi such that the two-dimensional space becomes stable. 
Next, one turns to the third cell of the matrix Eq. (13), 
puts P 2 = Pi, and implements the sign-changing algo¬ 
rithm complemented by augmentation of P 2 by unity, 
if needed, until the third dimension gets stable. The 
procedure being sequentially applied to all the cells of 
Eq. (131. finally yields a feedback matrix /3 stabilizing 
the required subspace of the dimension M. for the cho¬ 
sen task matrix a. 

For the case of N = 100, M = 10, in Table [2] we 
present sets of the negative real parts Sj of the 

eigenvalues Cj of the matrix 2/3 obtained after 10 cast¬ 
ings (to = 1... 10) for random matrices 2 and matrices 
/3 found as a result of implementation of the described 
algorithm. In Fig(2]we show the average eigenvalues and 
their mean absolute value deviation from the average 
calculated with the help of this algorithm. Both these 
quantities are exponentially decreasing with the eigen¬ 
value number. Note that it does not make sense to put 
the results of the calculations into Table |T| since the 
corresponding entries are all unities. 

This result means that the local control based on the 
change of the feedback sign for the most nimble vari¬ 
able (the AMN-rule) combined with a simple non-local 
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1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

{Si}, 

1.59 

1.57 

1.57 

0.00253 0.00253 

0.00161 

0.00161 

0.000144 

0.000144 

2.0 lO' 5 

2 

1.2 

0.219 0.219 

0.197 

0.197 

0.187 

0.0214 

0.0214 

1.69 lO" 3 

1.69 lO” 3 

3 

0.849 0.131 0.131 

0.0297 

0.0297 

0.0223 

0.0223 

1.51 10" 4 

2.3 10 -5 

2.3 lO" 5 

(Si} 4 

2.86 

1.28 

1.28 

0.317 

0.317 

0.0807 

0.0807 

5.75 10“ 3 

5.75 lO" 3 

2.26 lO” 3 

{A*} 5 

2.27 

0.37 

0.37 

0.0824 

0.0824 

2.46 10~ 3 

8.38 lO" 4 

8.38 lO" 4 

1.24 10 -4 

6.33 lO” 5 

{^} 6 

1.48 

0.614 0.614 

0.411 

0.13 

0.13 

0.109 

0.109 

0.1 

0.1 

{<^} 7 

3.54 

3.54 

1.21 

1.21 

0.208 

0.208 

0.0317 

0.0317 

2.87 lO" 3 

2.87 lO” 3 

{^} 8 

6.58 

0.751 

0.27 

0.27 

0.187 

0.187 

0.0363 

0.0363 

0.012 

0.012 


1.85 

1.85 

0.507 

0.507 

0.0478 

0.0478 

0.0244 

0.0244 

6.83 lO" 4 

6.83 lO” 4 

{*5j}io 

0.652 0.652 0.602 

0.602 

0.138 

0.138 

0.0734 

0.0734 

4.71 lO" 3 

4.71 lO” 3 


Table 2 Sets of the negative real parts of the eigenvalues of matrix a/3 obtained with the help of the feedback strength exploring 
algorithm for ten randomly chosen task matrices d. 


8 Re Si 



Re Si 


Fig. 2 The average sorted eigenvalues (solid line) and their 
average mean absolute value deviation (dotted line) obtained 
for the case N = 143, M = 11 with the help of the feedback 
strength-exploring algorithm. The average has been taken 
over 25 casts of the random matrix a. 


control solves the feedback search problem. This nonlo¬ 
cal control requires freezing the feedback parameters for 
the previous generations once a new generation starts to 
perform the AMN-rule-based local control. Moreover, it 
requires a learning procedure with adjusted attenuation 
of the generation couplings. Though the proposed algo¬ 
rithm solves the problem for any large dimension M, 
and thus explicitly demonstrates a possibility to con¬ 
struct a 100% successful algorithm, it may yield the 
coupling strengths so weak, that the obtained solution 
has no practical value. 


2.4 Susceptibility to Noise Reveals Feedback Structure 


We now address the question: What happens with the 
convergence towards the stabilized subspace ensured by 


the feedback matrix 

M 

f2ij = signet oo) ^2 Pi,kOtk,j (12) 

k =l 

in the presence of noise ? Again, as earlier, we empha¬ 
size that for physiological systems ’’noise” implies ex¬ 
ternal task-independent actions that can be of various 
origin, including both spontaneous changes in intrinsic 
neural variables and uncontrolled action of an external 
force field. Here signet —i oo) denotes the sign, which 
results from successful implementation of the control 
algorithm. We specify this question by asking: How far 
from the average positions Xi(t) satisfying the dynamic 
equations can the actual variables X, ; (f) = Xi(t) + 5xi(t) 
deviate in the presence of a time-dependent noise fi(t) 
? It can be immediately answered in the basi s o f the 
eigenvectors Xi of f2 it j where !&{t) = J2i x i(t)Xi- 
One considers the dynamic equations 

^Sxi(t) = Si8xi(t ) + ( 13 ) 

for the deviations 5xi(t), where S, denotes eigenvalues 
of f2. The solution 

t 

Sxi(t) = J e 5<(t_r) fi(r)dr, ( 14 ) 

o 

of this equations in the Fourier representation 
8xi(u) = 1 . ( 15 ) 

Oi — IUJ 

suggests that the spectral noise intensity |xi(w)| 2 of the 
coefficients Xi(t) is related to the spectral noise intensity 
\fi(u})\ 2 of fi(t) via the relation 

= — 1 . ,2 l/»MI 2 > ( 16 ) 

I Oi — IUJ I 

which corresponds to the susceptibility \Si — icu\ 2 . 
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In order to find the net mean square deviation yf (5xf) 

Jf duj |fei(w)| 2 /2tt of the variable X ) , one has to 
choose a model for the random perturbing forces or 
noise. In contrast to the Gaussian noise, typical of phys¬ 
ical multi-body systems, due to the central limit theo¬ 
rem of statistics, biological systems may, and usually do 
manifest other type of noisy behavior, since the under¬ 
lying processes may emerge from complicated nonlinear 
and non-random processes developing at the time scales 
much shorter as compared to the time scale of the pro¬ 
cess under consideration. The simplest model includes 
randomly chosen static forces /* acting on the relevant 
variables during a time interval T and, in each, next 
time intervals, these forces show random values. This 
yields 


I/*(<*>) | = 


sin 


ujT 


luT 


\h\\ 


and hence 
{5x1) = I fi 

= \h \ 2 


2 f did 


Sill 


ujT 


\Si — iu)\ X 

Re[(e* r -l)Sr 2 ] 

T Re Si 


(17) 


(18) 


Average square variances ) of the independent modes 
correspond to eigenvalues Ci of the covarience matrix, 
when one considers the problem in a basis, different 
from |d^j. 

For the case Si = 0 one should take into account a 
finiteness of the evolution time interval T and replace Si 
by i/T , thus obtaining the net mean squared deviation 
diffusively rising with time 


<**?> - r/T. 


(19) 


Noise analysis may serve as a powerful tool of re¬ 
vealing the eigenvectors of the matrix f2 r j of Eq. (12), 
which via the Eqs.( 16|17 ), also gives some ideas about 
magnitudes of the corresponding eigenvalues for a rea¬ 
sonably chosen model of the noise. One needs to cal¬ 
culate the noise covariance matrix GVy = (SxiSxj) from 
the experimentally observed deviations of the running 
measured values from their ave rage values in the sta¬ 
tionary regime. Eigenvectors Xi of this matrix should 
coincide with those of 17 = G(t —► oo)/3a. Moreover, 
this method can also be employed for determining the 
eigenvalues of the time-dependent values of the feed¬ 
back matrix Q(t) = G(t)/3a , provided the time intervals 
where this object changes significantly are long enough 
and exceed a typical correlation time of the noise. 

Note that for the case N > M, the non-zero eigen¬ 
values of the N x N matrix i? = Gf3a coincide with M 


=eigenvalues of the MxM matrix cD = aG/3. The remain¬ 
ing N — M eigenvalues are zeros, unless an additional 
requirement is imposed on the matrix fl. Spread of the 
variations 5xf in the directions of the corresponding 
eigenvectors are expected to be large and be governed 
by the spectral properties of the noise. In other words, 
given a task matrix and a feedback matrix, which has 
been found with the help of the stabilization search 
algorithm, there exist two subspaces in the space of 
the variables X{, the so-called uncontrolled manifold 
(UCM) [57] and its orthogonal complement (ORT) also 
sometimes addressed as the null-space and the range- 
space, respectively. ORT is the task-specific subspace 
expected to show high stability, which implies that vari¬ 
ance in ORT is expected to be small. In the UCM, vari¬ 
ance is generally expected to be large and increasing 
with time, unless there are other factors, outside the 
explicit task formulation, that keep it within a certain 
range. In fact, after a person becomes an expert in a 
motor task, repetitive attempts at the task lead to vari¬ 
able combinations of elemental variables that all lead to 
successful task performance within a permissible error 
margin (0; reviewed in mi). However, not all pos¬ 
sible combinations of the involved elemental variables 
are used across repetitive attempts. Such self-imposed 
additional constraints may be addressed as ’perfection¬ 
ism’; they may reflect optimization with respect to a 
cost function (e.g., [36]). Figure [3] illustrates the task 
Xi+ X 2 = C for different values of C. While all points 
on the slanted dashed lines correspond to perfect task 
performance, actual behavior shows much more con¬ 
strained clouds of solutions that show larger deviations 
along the solution space (the UCM) compared to devi¬ 
ations orthogonal to that space (that lead to errors in 
performance). Note that if the task is learned for a par¬ 
ticular value of C, the solutions show robustness: Sim¬ 
ilar relative locations and shapes of the solution clouds 
for other values of C. 

More formally, assuming a feedback matrix /3 sta¬ 
bilizing a M --dimensional subspace is found, one may 
think of imposing complementary constraint, that can 
either be in the form of explicit equations or follow from 
a requirement of minimization of a cost function. In 
such a setting, one speaks about stabilization of a sub¬ 
space of a dimension M + exceeding the dimensionality 
M of the initial task subspace. Such perfectionism may 
be viewed as a secondary tasks decreasing variance in 
some directions of the UCM. 

We would like to emphasize once again that there 
exists three types of specific bases that can be employed 
for the dynamics description: (1)A laboratory basis of 
measured variables, (2) A basis of modes resulting from 
both body- and task-imposed constraints, while each 
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Fig. 3 Consider a task of producing a constant sum of two 
variables, e.g., pressing with two fingers and producing a con¬ 
stant total force level. The dashed lines show solution spaces 
(uncontrolled manifolds, UCMs) for three different force lev¬ 
els, Ci, C2, and C3. Across repetitive trials, clouds of data 
points form ellipses elongated along the corresponding UCM. 
This shape reflect lower stability along the UCM as compared 
to the orthogonal direction relevant to the task-imposed con¬ 
straint. Note that the three data clouds are centered not ran¬ 
domly along the UCMs but reflect a certain preferred sharing 
of the task between the two effectors. This preference may 
reflect an optimization principle. 


mode has independent noise statistics, and (3) A task- 
specific basis of variables along which the control is lo¬ 
cal. The last basis is a conceptually new feature afford¬ 
ing a simple structure of the control algorithm. 

The relation between the second and the third bases 
may change during the process of stability search: Ap¬ 
plication of the local control G leads to a change of 
Q = aG(3, that is, the way the task requirement a is 
mapped onto the feedback action. Thereby it may af¬ 
fect both the basis of the noise-decorrelated modes and 
the magnitude of the corresponding modes susceptibil¬ 
ities. There exists such an extreme, where the mode 
susceptibilities |«S* — iui\ 2 vary without changing the 
corresponding eigenvectors Xi of 17 - the linear com¬ 
binations ~?t\ and X j of the laboratory basis variables 
remain statistically independent, and only the variance 
(Sxf) of one combination starts to exceed that of the 
other (Sxj). One can call this regime ’’mode crossing”, 
by the analogy to the phenomenon of ’’term crossing” 
in Quantum Mechanics [T]. In the general case where, 
along with a change of the susceptibilities, the local 
changes equally result in the emergence of an apprecia¬ 
ble covariance between the modes ( SxtSxj ) ^ 0 , one en¬ 
counters the phenomenon of so-called ” avoid crossing”, 


where formerly the larger mode susceptibility, though 
approaching in magnitude the other_one, remains al¬ 
ways larger, while both eigenvectors X j and X :j corre¬ 
sponding to these modes rotate. 

2.5 Algorithm for Hierarchical Feedback 

Hierarchical control considered in this section gives an 
example of stabilization search algorithm different from 
that considered earlier. We demonstrate that such an 
algorithm requires only local control, whereas the role 
of nonlocal control can be played by another randomly 
chosen linear feedback matrix once the current one does 
not yield stabilization. The idea of hierarchical control 
of the human body is very old. Arguably, the first hi¬ 
erarchical system considered the brain and the spinal 
cord. A comprehensive scheme of control with referent 
configurations has been suggested recently built on a 
hierarchical principle, starting with referent values for 
a few task-specific, salient variables, and resulting in 
referent length values for numerous involved muscles 
|l 6 ) . An example of hierarchical control is the command 
structure of an army, where a general gives orders to 
privates through the commanders of regiments, battal¬ 
ions, and platoons. We explore efficiency of hierarchical 
control in this section, along with another modification 
- the most nimble x * is no longer punished by a step 
change of the sign of its contribution but experiences a 
smooth change of the participation in the feedback as 
the cosine of the corresponding functional z t . 

2.5.1 Intrinsic instability of the hierarchical control 

In mathematical terms, the time derivative of a vector 
n of variables at ro-th step depends on the variables 
"af n _iat the previous step, while the spatial dimension 
N n of each next step varies from step to step. At each 
step local control may be implemented. For example, 
the corresponding set of the local control equations for 
a three-level hierarchy has the form 

= 3^3 ( 20 ) 

5 = ( e 0 „ 

N„( S H 

| a » i =(6a) u (e* 2 ) i , 

based on the diagonal matrix elements 
G n ) = cos [z n i (f)] of the local control operators G n . 

If at each level n of the hierarchy, we also include di¬ 
agonal damping matrices 7 n of the dimension given by 
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the number N n of the variables at this level, the set of 


equations (20) can be seen as a single matrix equation 


71 ) 

^t' 

\ 




72 ) 

^2 

= 




%) 

^3, 

I 





1 Gi 

0 

0\ 

(o 

0 Aa 


0 

g 2 

0 

B 

0 0 


V 0 

0 

g 3 J 

1° 

C 0 



( 21 ) 


Alternatively, the local control operation may be ap¬ 
plied not to susceptibility of a given variable to external 
factors, but to efficiency with which this variable acts 
on the variables of the next generation. In such a case, 


two operators on the left hand side of Eq.(21) have to 
be interchanged 

'(&- 7i) Ji'. 

( i - 72)^2 = 


4 - 73) ^ 


k dt 


/ 0 0 as n 


BOO 


v 0 C 0 , 

V 



( 22 ) 


Though offering a simple way to address many vari¬ 
ables at once, hierarchical control may add instability. 
Therefore, the intermediate steps have to be damped 
in order to avoid such an instability in contrast to one- 


step control of Eq.(10). In particular, note that for zero 


damping 7 ,; = 0, even the simplest two-level control be¬ 
comes unstable, and this is always the case for a higher 
number / of the control levels. The root of this insta¬ 
bility is rather simple, and can be illustrated with an 
example of a three-level control scheme, with just one 


variable at each level, when Eq. (21) takes the form 
£ t xi = GiAaxz 


d 

if 

jr t X 3 = G 3 CX 2 


2 = G 2 Bx 1 


(23) 


and hence the roots of the characteristic polynomial 
given by /-th roots of the eigenvalues of the matrix 
Gi... G 3 CG 2 BGiAa are also uniformly distributed on 
the corresponding circles in the complex plane, with 
radii given by the eigenvalues moduli. 

We thus come to the conclusion that damping is in¬ 
dispensable for stable hierarchical control. At the same 
time, all 7 * > 0 imply a trivial asymptotic situation 
of completely damped motion at all levels for all tasks. 
Therefore, in order to have a reasonable model, we have 
to assume vanishing of the damping rates at only one hi¬ 
erarchical level of control. Presumably, this level should 
have the maximum dimensionality N n . 


2.5.2 Action of the feedback loops 

We note that stability can be to a certain extend im¬ 
proved by introducing intermediate feedback loops to 
the net feedback loop as shown in Fig|4j 



corresponding to the characteristic equation 

A 3 = G 3 CG 2 BG\Aa. (24) 


It is evident that, whatever a non-zero complex number 
standing on the right hand side of this equation is, the 
phase factors of the three roots of the cubic equation 
are equally distributed on the unit circle in the com¬ 
plex plane such that at least one of them has a positive 
real part. The same structure of the root distribution 
persists in a higher-dimensional case with /-level con¬ 
trol, since the characteristic equation in this case has 
the form 


A 1 -Gi... G 3 CG 2 BGiAa 


= 0 , 


(25) 


Fig. 4 A scheme of hierarchical control of the hand within 
the idea of control with referent configurations (RCs) of the 
body [TO]. At the top level, a low-dimensional set of referent 
values for salient, task-specific variables is reflected in the RC. 
A sequence of few-to-many transformations results in higher¬ 
dimensional RCs at the digit level and muscle level. Local 
feedback loops ensure stability with respect to the variables 
specified by the input. The global feedback loop ensures that 
the actual body configuration moves towards one of the solu¬ 
tions compatible with the task RC. At each level, inputs to 
a neuronal pool (N1,N2, and N 3) are combined with affer¬ 
ent feedback (AFF) to produce the output (efferent signals, 
EFF). At the lowest level, elements are alpha-motoneurons 
and their referent coordinates correspond to the thresholds 
of the tonic stretch reflex (lambda). Modified by permission 
from pL6| . 


Det 
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For the example of aforementioned three-level hier¬ 
archical control, the Laplace-transformed equations of 
motion for the case of feedback loop equipped with sec¬ 
ondary feedback loops read 

Alz^i = G\Aont 3 + d. 2^2 
A ^ 2 = G 2 Bltf i + A 3 lt 3 

XHcf 3 = G 3 C^2 (26) 

where the secondary feedback loops of the second and 
the third hierarchical levels are given by the matrices 
A 2 and A 3 , respectively. The latter can be also consid¬ 
ered as dependent on local variables z. For a given set 
of variables z, the instantaneous eigenvalues A can be 


or the corresponding variable velocity, 


calculated as the roots of an analog of Eq.(25), which 
reads 


Det 


1 


1 


A - G 3 C _ _ G 2 B _ _ GiAa 

A — G 2 BA 3 A — G\AaA 2 


= 0 . 


Note that in the basis where the matrices G 2 BA 3 
and G\AaA 2 are diagonal, they can be interpreted as 
damping matrices 71 and 72 , respectively, although in 
this case these matrices are dependent on parameters z, 
which, generally speaking, in this representation cannot 
be considered as local control. Still, this brings about an 
idea of the self-adjusting damping dependent on time 
via the parameters Zj(f). 

2.5.3 Stabilization by self-adjusting damping 

Adjusting gain in the feedback loops controls stability 
at each hierarchical level. This effect can be modeled 
when the damping parameters 7 i are taken depending 
on the local parameters Zj, increasing, for instance, as 


7 i = Bi+ iMZi 


(27) 


, for the variables at each damped level. The parameters 
7 j are positive for all but one hierarchical levels, where 
they may remain zero, in order to avoid the trivial case 
of complete damping of all variables. An analog of the 


characteristic equation (25) for this case gets simplified 
and reads 


Det 


A - 73 - G 3 C 


1 


G 2 B- 


1 


A - 72 A - 71 


-Gi Aa 


= 0 , 


where one of the three diagonal matrices 71 , 72 or 73 is 
assumed to be zero. 

In turn, for the local control over the damping rate, 
the dependence Zi(t ) can be given by a differential equa¬ 
tion relating the positive rate of increase f^Zi with an 
even power / of either the corresponding variable, 


dt 


Zi = e (xi) f , 


(28) 


dt Zi = £ 


dxi 

dt 


f 


(29) 


where e < 1 is a numerical parameter, which may de¬ 
pend on the hierarchy level number and the variable 
number. 

In the following Table [3] we present the success rates 
calculated for series of 100 random casts of the rectan¬ 
gular feedback matrices A, B, and C corresponding to 
the dimensions N 3 , N 2 , and N 3 of the spaces of x\, x 2 , 
and x 3 , respectively, and generated independently each 
time for a random cast of the N 3 x M task matrix a. 
The identical damping parameters Tj = 0.1 and /x, = 1 
have been chosen. Moreover, in order to avoid the triv¬ 
ial situation yielding 100 % success rate, where all vari¬ 
ables get strongly damped and hence become vanish¬ 
ing in the course of time, a fourth hierarchical level 
corresponding to the vector variables x 0 defined by 
the differential equation without damping, 1 lif 0 = l/j 


has been included into Eq. (231. The non-zero diagonal 


matrix elements of the local feedback control matrices 
Gi, 2,3 have been chosen here as Ga = cos( 0 . 1 z*), and 
the power index / and the numerical parameter e in 


Eq. (29) have been chosen equal to 4, and 0.1, respec¬ 


tively for all variables. 

A typical time dependence l/ (t) for the successfully 
stabilized subspace is shown in Fig{5] (the top plot). In 
accordance with the smoothness of the local feedback 
control dependence (cos(O.lzj) instead of sign(107r — 
Zi)), the time dependence does not manifest a cusp-like 
changes of the derivatives, typical of those seen in Figjl] 
For unsuccessful control, the dependence looks like the 
one depicted in Fig(5] (the bottom plot). 

A rather high success rate of the hierarchical local 
control over the dynamics for random choice of the feed¬ 
back matrices and self-adjusting damping suggests a 
strategy, which can replace the non-local control. Once 
the current random cast of the feedback matrices does 
not yield stabilization during a run, one may achieve 
it by taking another random cast in the course of the 
same run, and in case of failure, repeat the attempt 
again and again. In particular, the probability 0.3 im¬ 
plies that just a few (4 — 5) such casts during a run are 
required in order to stabilize a 5-dinrensional subspace. 
Changing the entire linear feedback matrices when the 
local control turns out to be inefficient may be consid¬ 
ered as an action replacing nonlocal control discussed 
in the previous section. 

Note, that the dimensionality at each step of the hi¬ 
erarchy does not need to be larger than the dimension¬ 
ality of the previous step. We illustrate this by an ex¬ 
ample of control over a two-dimensional space y exerted 
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M 

1 2 3 4 5 6 

N 1 /N 2 /N 3 

5/50/200 5/50/200 5/50/200 5/50/200 5/50/200 6/50/200 

R 

1 1 0.69 0.43 0.3 0.15 


Table 3 The stabilization success rate of a Af-dimensional subspace for hierarchical control with Ai/^/lVa-dimensional levels. 


2 -object dimension 16 -brain dimension 


4 -finger dimension 








Fig. 6 Examples of successful (lower panels) and unsuccessful (upper panels) search for the new equilibrium control. We 
show coordinate variables as functions of time in arbitrary units. First column: deviations of two controlled coordinates from 
their prescribed value - zero. The second column: dynamics of the hypothetical elemental variables tend to some asymptotic 
values for the successful control. Third column: dynamics of motor variables. One sees the dynamics of transition from the old 
equilibrium position to the new one. 


by a three-step hierarchical feedback with locally chang¬ 
ing parameters. The first-step variables x \ belong to 
a space of the dimension N\ = 4; here the damping 
increases according to Eqs.( 27|29 | with e = 0.2, / = 2, 
and Hi = 0.1. Next step of the hierarchy comprises the 
variables it 2 that belong to a space of the dimension 
N -2 = 16 where they experience no damping. The last 
step of the hierarchic feedback comprises variables it 3 
in a four-dimensional space, N :i = 4, with strong con¬ 
stant damping T) = 5 at the third level. The success 
rate for these parameters was R ~ 0.8. Examples of un¬ 
successful and successful search of stability are shown 
in Fig[6]for the components of the vectors y, it 2 , and 
X 3 . One sees that for successful control the components 
of the high-dimensional vectors tend to asymptote with 
time, while for unsuccessful control, they keep chang¬ 
ing. 


3 Experiments with Human Hand 

An experiment was performed to illustrate one of the 
central ideas of the suggested scheme and check some 
of its predictions. We used the task of accurate force 
and moment-of-force production by the four fingers of 
the dominant hand. The forces produced by four fin¬ 
gers, the index, the middle, the ring, and the little, 


were controlling the position of a point on the computer 
screen. Two types of perturbations were used. First, we 
modified the visual feedback leading to changes in the 
mapping between the finger forces and the two task 
variables. Second, we used mechanical perturbations 
applied to a finger that led to actual changes in the 
two task variables. The experiments were approved by 
the Office for Research Protections at the Pennsylvania 
State University. 


3.1 Visual Feedback Perturbation 

The first experiment was as follows. Four 6-axis force/moment 
sensors (Nano-17, ATI Industrial Automation, USA) 
mounted on the table were used to measure normal 
forces produced by the tips of the index (/), middle 
(M), ring ( R ), and little (L) fingers. To increase fric¬ 
tion between the digits and the sensors, 320-grit sand¬ 
paper (SandBlaster,3M, USA) was placed on the con¬ 
tact surfaces of the sensors. The centers of the sensors 
were evenly spaced at 30 mm. The output analog sig¬ 
nals from the sensors where digitized with the 16-bit 
data acquisition card (PCI-6225, National Instrument, 
Austin, TX, USA) at 100 Hz. A Lab VIEW program 
(LabVIEW 2011, National Instrument, USA) was used 
to provide visual feedback and store the data on the 
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Fig. 5 Time dependencies for the successfully controlled five¬ 
dimensional subspace in the hierarchical cascade setting with 
N 1 /N 2 /N 3 = 5/50/200 (upper). Unsuccessful control (lower). 


computer (Dell Inc., USA). Offline processing and anal¬ 
ysis was done in Matlab (MathWorks, USA). 

At the beginning, the four-dimensional space of the 
finger forces was mapped onto the two-dimensional space 
of the screen according to a very natural rule - the 
vertical coordinate was proportional to the net resul¬ 
tant finger force, while the horizontal displacement was 
given by the net moment of the finger forces computed 
with respect to a horizontal axis passing in the anterio¬ 
posterior direction in-between the R and M finger sen¬ 
sors. The subject was first requested to place the cursor 
into a position in the middle of the screen and, second, 
to keep it in this position at all times. Once this task was 
accomplished, the law according to which the deviation 
of the finger forces from the steady-state finger force 
values are mapped to the deviation of the point from the 
center of the screen was changed without the subject’s 
knowledge. The new law relating the four-dimensional 
space of the force deviations with the two-dimensional 
space of the cursor deviations from the center point was 
given by a new, randomly chosen, 2x4 Jacobian ma¬ 
trix. This law drew the system to an unstable regime. 
In order to keep the point at the center of the screen, 
the subject had to find a new feedback rule stabilizing 
this position. Trajectory of the forces exerted by the fin¬ 
gers in the course of this search for new stability were 
recorded. 


3.1.1 Results - manifestations of the feedback changes 

In FigjT] we present an example of the time profiles of 
the finger forces (upper panel), along with the coordi¬ 
nates of the point on the screen (lower panel), which 
were recorded during a successful trial. Overall stabi¬ 
lization success rate was R ~ 0.55. 



Fig. 7 Finger forces in N as functions of time in 10 _2 s (upper 
panel) and the coordinates of the point position on the screen 
(lower panel). The law relating the forces to the position of 
the point was changed at t = 7s, and the values of the finger 
forces at that time were chosen as zero N. 


There is a qualitative similarity between the depen¬ 
dencies depicted in FigjT] and the corresponding calcu¬ 
lated dependencies illustrated in Fig[6| 

The sharp spikes and jumps on the experimental 
curves correspond to changes of the trends; these events 
show approximately the moments of time when a cor¬ 
recting action was taken similarly to the simulated curve 
in FigjT]. Though spikes in the simulated curves may be 
unambiguously associated with the correcting actions 
of the feedback sign change, within the experimental 
finger force curves analysis there is no formal rules al¬ 
lowing to identify such moments, and one can speak 
only about an intuitive similarity between the depen¬ 
dencies. Moreover, since we are exploring the regime 
of searching for an equilibrium in new, formerly un¬ 
known, conditions, we cannot invoke the powerful tool 
of statistical analysis comparing different trials, since 
each new trial corresponds to new initial and task con¬ 
ditions, whereas multiple sequential repetitions of the 
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same experimental setting for the same subject would 
likely involve processes of learning and adaptation that 
are beyond the scope of questions addressed here. Still 
certain information about the moments of changes of 
the feedback matrix can be extracted from the analy¬ 
sis of noises, since, as discussed in Sect |2.4[ the princi¬ 
ple axes of the tensor susceptibility to noise (so-called 
principle components) coincide with eigenvectors of the 
dynamic matrix 17, and when directions of these axes 
change, the matrix 17 changes as well, thus implying a 
change of the feedback matrix. Moreover, the larger the 
eigenvalue Ci of 17 is, the smaller is the susceptibility 
of the corresponding direction to noise, in accordance 
with Eq.( 16). The noise analysis is discussed in the next 
section. 


they are decreasing exponentially with each next eigen¬ 
value being scaled, on average, by a factor, as it was 
the case for the hierarchy of nimbleness shown in Fig(2] 
This has a consequence important for the noise analysis: 
The directions that belong to the two-dimensional sub¬ 
space relevant to feedback are less susceptible to noise, 
while the noisy directions correspond to the null-space 
and do not contribute significantly to the feedback. 

The dependence Xi (n) of the forces produced by the 
I, M, R, and L fingers was captured at the sequential 
time points t = nr separated by time intervals of r = 
10 -2 sec. The covariance matrix was extracted from 
these experimental data in several steps. First, for each 
Xi(n) an average time dependence Xi(n) = Xi(t) has 
been calculated numerically as 


3.1.2 Analysis of the principal components of noise 
covariance reveals feedback changes 


x i{t) = 

n 


Xj{n) 

y -\/tt 


exp 


(t ~ nf 

Y 2 


(31) 


Human fingers are not independent force generators: 
When a person tries to press with one finger, other 
fingers of the hand show unintentional force increase 
mm- This phenomenon has been addressed as en¬ 
slaving or lack of finger individuation Patterns 

of enslaving are person-specific and relatively robust; 
changes in these patterns have been reported with spe¬ 
cialized practice [3U ?]■ These patterns may be described 
as eigenvectors of enslaving E,i = {fj}& where i = 
{I — index; M — middle, R — ring, andL — little} stands 
for an instructed finger. Directions of E,i may be 
viewed as preferred directions in the space of finger 
forces when the person is trying to press with individual 
fingers They are related as 

1 f E,i = uli (30) 

to the forces X, of the individual fingers i by an or¬ 
thogonal matrix U representing rotation in the four¬ 
dimensional space. We assume, that these very direc¬ 
tions may change during the search of stability and they 
manifest statistically independent fluctuations thus be¬ 
ing the eigenvectors A, of the deviation covariance ma¬ 
trix. 

Since in our experiment the task was formulated in 
a two-dimensional space, and hence the rank of the dy¬ 
namic problem presumably equals two, there should be 
only a two-dimensional sub-space in the space of the 
finger forces that governs the dynamics of the point on 
the screen. This implies that only two out of four eigen¬ 
values of the matrix 17 differ from zero, and the other 
two vanish, unless an external to the task requirement 
associated with additional constraints is imposed upon 
the system, as mentioned at the end of Sect |2.4| In the 
latter case, all the eigenvalues do differ from zero, but 


where r is chosen as a time unit, while the averaging is 
performed over a time window ~ Y r with a width Y. 
Next, the covariance of the finger forces 


CiAt) = T, 

n 


Sxi(n)Sxj(n) 

Yyfn 


exp 


{t~nf 

Y 2 


(32) 


where 


Sxi(n) = Xi(n ) - Xi(n), 


was found numerically with the same Gaussian width 

Y. 

Being a real and symmetric matrix, GY ? can be set 
in a diagonal form by an orthogonal transformation 
given by a rotation matrix Uij (t) and its inverse matrix 
UrHt), such that 

= y 'UiMWkWkAt). (33) 

k 


The eigenvalues Ck (t) provide us with the principle com¬ 
ponents of the noise in the orthogonal directions of 
statistically-independent modes, while the matrix Ui t k(t ), 
which in the laboratory basis can be viewed as a row of 
the column eigenvectors Xi{t ), relates these modes to 
the individual finger variables. All these quantities were 
found numerically from the data obtained for C, : j(t). 
Note that thus obtained orthogonal matrix U i: k{t) ex¬ 
periences a time evolution corresponding to rotation in 
the 4-dimensional space of the finger forces, while the 
angular velocity of this rotation can be found as the 
eigenvalues of the left logarithmic derivative of Ui^it) 
defined as 


R{t) = - d ^U~\ 


(34) 
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The eigenvalues Ri of R(t) are real and have pairwise 
opposite signs, such that only two real numbers charac¬ 
terize rotation in the four-dimensional space. We calcu¬ 
late these quantities replacing the derivative in Eq. (341 
with the finite difference between two neighboring in¬ 
teger time points. In order to get rid of the so-called 
shot noise, which is an error-inducing influence of such 
a replacement, the calculation must be followed by av¬ 
eraging over a time interval shorter than Y. 


Results of such processing of the experimental data 
presented in Fig(7]are shown in Fig|8] 

One can identify eleven rotations of the covariance 
matrix basis Xi{t) presumably associated with changes 
of the feedback matrix. Note that the highest rotation 
velocity does not necessarily produce a strong effect on 
the finger forces, since the relevant quantity rather cor¬ 
responds to the spike area, representing the rotation 
angle. The situation has much in common with dynam¬ 
ics of so-called adiabatic and diabatic molecular term 
crossings, well-known in Quantum Mechanics[lj. The 
maximum rotation velocity corresponds to the time mo¬ 
ments when two or several eigenvalues of the matrix 
have a tendency to coincide, thus getting rid of the dif¬ 
ference between the large noise typical of the null-space 
and small noise typical of the orthogonal subspace in a 
stationary regime. 


3.2 Experiments with Mechanical Perturbations 


The results of the first experiment have demonstrated 
qualitative consistency of the model and the experi¬ 
ment. Still the main assumption of the suggested stable 
control search algorithm, the principle “act on the most 
nimble one” (AMN) requires additional arguments. The 
second experiment was designed to test this very prin¬ 
ciple. Mechanical perturbations (lifting and lowering 
a finger) were applied during the performance of an 
accurate multi-finger steady-state task. According to 
our scheme, quick reactions to these perturbations are 
based on the AMN rule. We checked this prediction 
by comparing the directions of changes in the finger 
force space produced by unexpected perturbations of 
the steady-state force patterns (described as a vector 
7 p,i ) with the first identifiable correction produced by 
the subjects (a correction vector, c,i)- We expected 
the angle between / c,i and to be small, smaller 
than the angle between j c,i and E,i, defined by 


Eq. (30). 



Fig. 8 The angular velocity R\ (t) of rotation of the basis of 
the covariance matrix eigenvectors as function of time (bot¬ 
tom, solid line, arbitrary units). Logarithms of the eigenval¬ 
ues Cfc(t) of the covariance matrix (four bold curves above 
the velocity). The average of the covariance matrix was com¬ 
puted over Y = 100 sequential time points with the interval 
of 10 -2 sec between the points. The switching of the feed¬ 
back eventually took place at the domains of ’’avoid cross¬ 
ings” discussed at the end of Sect |2.4| marked here by arrows, 
and corresponding to the maxima of the angular velocity. On 
the top of the plots, the finger forces and the cursor coordi¬ 
nates, corresponding to the dependencies in Fig[7] averaged 
over the same time intervals Y = 100, are shown for com¬ 
parison. The ’’avoid crossing” of the three eigenvalues occur 
for the 4-th, 7-th, and 8-th intersections. The strongest con¬ 
tribution comes from the 5-th crossing, which presumably is 
relevant to changing of the direction in the two-dimensional 
orthogonal subspace. 


3.2.1 Methods - the ’’inverse piano” setup 

Eight young, healthy subjects took part in the exper¬ 
iment (four males). They were right-handed, had no 
specialized hand training (such as playing musical in¬ 
struments) and no injury to the hand. 

An “inverse piano” apparatus was used to record 
finger forces and produce perturbations. The appara¬ 
tus has four force sensors placed on posts powered by 
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linear motors, which could induce motion of the sen¬ 
sors along their vertical axes (for details see Martin et 
al. 2011). Force data were collected using PCB model 
208C01 single-axis piezoelectric force transducers (PCB 
Piezotronics, Depew, NY). The signals from the trans¬ 
ducers were sent to individual PCB 484B11 signal con¬ 
ditioners - one conditioner per sensor - and then dig¬ 
itized at 1 kHz using a 16-bit National Instruments 
PCI-6052E analog-to-digital card (National Instruments 
Corp., Austin, TX). Each sensor was mounted on a Lin- 
mot PS01-23x80 linear actuator (Linrnot, Spreitenbach, 
Switzerland). Each actuator could be moved indepen¬ 
dently of the others by means of a Linrnot E400-AT 
four-channel servo drive. Data collection, visual feed¬ 
back to the subject, and actuator control were all man¬ 
aged using a single program running in a National In¬ 
struments LabView environment. Visual feedback was 
provided by means of a 19” monitor placed 0.8 m from 
the subject. The feedback cursor (a white dot) showed 
on the monitor represented total finger force along the 
vertical axis and total moment of force in a frontal 
plane computed with respect to a horizontal axis pass¬ 
ing through the midpoint between the M and R fingers 
along the horizontal axis. Pronation efforts led to left¬ 
ward deviation of the feedback cursor. An initial target 
was placed on the screen (a white circle) corresponding 
to the total force of 10 N and zero total moment. 

The experiment involved two parts: Voluntary force- 
pulse production and reacting to unexpected perturba¬ 
tions (see later). Prior to each trial, the subjects were 
asked to place the fingers on the centers of the sensors 
and relax; sensor readings were set to zero during this 
time interval. As a result, the sensors measured only ac¬ 
tive pressing forces. Then, a verbal command was given 
to the subject and data acquisition started. The sub¬ 
ject was given 2 s to place the cursor over the initial 
target. During the force-pulse trials, the subjects were 
asked to produce a force pulse from the initial target in 
less than 1 s by an instructed finger (Figure [9]4). Each 
finger performed three pulse trials in a random order. 
In perturbation trials, one of the sensors unexpectedly 
moved up by 1 cm over 0.5 s. This led to an increase 
in total force (Figure [9^3) , while changes in the mo¬ 
ment of force depended on what finger was perturbed. 
The subject was instructed to return to the target po¬ 
sition as quickly as possible. Each finger was perturbed 
once per trial, with 10 — s rest periods between each 
of the three repetitions. Perturbation conditions were 
block-randomized between fingers with 1-min rest pe¬ 
riods between blocks. 

From the pulse trials we manually identified the on¬ 
sets of the force decrease phase (see Figure [9jA) . From 
the perturbation trials we manually identified the on- 




Fig. 9 Typical time profiles of the force-pulse trial (A, top 
panel) and perturbation trial (B, bottom panel) performed 
by a representative subject. The total force profile is shown 
with the solid trace, and individual finger forces are shown 
with different dashed traces. The intervals used to compute 
the eigenvectors in the force space are shown with the arrows. 


sets of two time intervals: one corresponding to the 
perturbation-induced force change and the other cor¬ 
responding to the earliest corrective action by the sub¬ 
jects. Each time interval contained 200 ms of the four¬ 
dimensional finger force (/, M, L, and R) data. Further, 
for each time interval, principal component analysis 
(PCA) based on the co-variation matrix was used to 
compute the first eigenvector in the four-dimensional 
finger force space, accounting for most variance across 
the time samples, for each subject and each trial sepa¬ 
rately. 

Therefore, for each perturbation trial we obtained 
two eigenvectors in the finger force space. We will refer 
to these vectors as "7* p^ (force vector during the per¬ 
turbation applied to the i-th finger), and c,i (force 
vector during the earliest correction in trials with per¬ 
turbations applied to the i-th finger). For each finger, 
from the three force-pulses trials we computed an av¬ 
erage vector / e,i (force vector during the downward 
phase of force change in the pulse task by the i-th fin- 
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ger). Note that this vector reflected the unintentional 
force production by non-task fingers of the hand (en¬ 
slaving) . 

Finally, we computed the angles apE (averaged across 
repetitions) between py and / E,i and angles ape be¬ 
tween ~J?p,i and c,i for each finger and each subject 

separately. The Harrison-Kanji test, which is an analog 
of two-factor ANOVA for circular data, was used with 
FINGER (4 levels: /, M, R, L) and ANGLE (2 levels, 
apE and ape ) as factors. All data analysis were per¬ 
formed in Matlab (MathWorks, Inc.) software. 

3.2.2 The results - acting along the most nimble 
direction 

During the force-pulse trials, forces of all four fingers 
changed in parallel (Figure There was a larger 
change in the force produced by the instructed finger 
and smaller changes in the other finger forces. These 
patterns are typical of enslaving reported in earlier stud¬ 
ies [44ll6] . PCA applied to the finger force changes pro¬ 
duced similar results over the phase of force increase 
and the phase of force drop. The first PC accounted for 
over 95% of the total variance in the finger force space 
in all subjects and for each finger as the instructed fin¬ 
ger. The loading factors at individual finger forces were 
of the same sign. 

In the perturbation conditions, lifting a finger’s force 
sensor produced a complex pattern of changes in the 
forces produced by all four fingers (similar to the re¬ 
sults described in mm- Typically, the force of the 
perturbed finger increased, while the forces produced 
by the three other fingers dropped (Figure The 
total force increased. The first PC accounted for over 
95% of the total variance in the finger force space in 
all subjects and for each finger as the perturbed finger. 
The loading factors at different fingers were of differ¬ 
ent signs; most commonly, the perturbed finger loading 
was of a different sign as compared to the loading of 
the three other fingers. 

Overall, the angle between the vectors of perturba¬ 
tion and correction ( ape ) was consistently lower than 
the one between the vectors of perturbation and volun¬ 
tary force drop (olpe)- These results are illustrated in 
Figur^Toj which shows averaged (across subjects) values 
of the two angles with standard error bars. The gross 
average of ape was 15.9±6.6°, while it was 25.9± 10.9° 
for apE- The Harrison-Kanji test confirmed the main 
effect of angle (F[- t i56 ] = 19.17,p < 0.0001) without 
other effects. 

Overall, these results confirm one of the predictions 
of the AMN-rule. Indeed, the first reactions to per¬ 
turbations in the four-dimensional finger force space 



Finger 

Fig. 10 The angles between the force vector produced by a 
quick perturbation applied to a finger and the force vector 
during the downward phase of the force-pulse trial by the 
same finger (app) and between the first vector and the vector 
of the corrective action ( ape )• Averaged across subjects data 
are shown with standard error bars. Note that for each finger 
as the target finger app < ape- 


showed relatively small angles with the vector reflect¬ 
ing the effects of the perturbation on finger forces. This 
angle was significantly smaller as compared to typical 
subject-specific finger force vectors produced when the 
subjects tried to increase or decrease force with one fin¬ 
ger (the perturbed finger) only. 


4 Discussion: The Stability Search Algorithm 
Hypotheses within the Context of Motor 
Control Hypotheses 

The most important axiom forming the foundation of 
our approach is that we assume the existence of task- 
specific coordinate systems organized to allow effective 
local control. A particular implementation of local con¬ 
trol has been addressed as the “act on the most nim¬ 
ble” (AMN) rule. We have shown that this method can 
solve problems better than control with random matri¬ 
ces but loses efficacy with an increase in the task dimen¬ 
sionality, not so much with the system dimensionality. 
This problem can be overcome by using feedback with 
adjustable gain, but in this case the system slows down 
dramatically. 

Further, we considered a number of additional rules 
that improve the outcome. One of them is: If local con¬ 
trol does not work, change the coordinate system. More 
specific rules that all improve probability of reaching 
stability include the following. Deal with one dimension 
at a time and do not return back to any of the previ¬ 
ously involved dimensions. Organize elemental variables 
into generations (groups assembled by links to a specific 
task variable, y) by their nimbleness. Allowing bi-local 
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control improves the performance even more. Bi-local 
control is an example of implementing the general rule 
of minimizing non-local actions. 

4.1 Systems of Coordinates in Motor Control 

One of the important features of the suggested scheme 
is the identification of three systems of coordinates that 
can be used to describe processes associated with the 
neural control of a movement. Most commonly, move¬ 
ment studies operate with variables directly measured 
by the available systems, for example those that mea¬ 
sure kinematic, kinetic, or electromyographic variables. 
Some of these variables describe overall performance, 
for example fingertip coordinates during pointing. Other 
variables reflect processes in elements that contribute to 
the task-related performance (e.g., joint rotations, digit 
forces and moments, muscle activation, etc.). Some vari¬ 
ables may not be directly measured but computed based 
on other variables and known (or assumed) parameters 
of the system (for example, joint torques and muscle 
forces). 

Using measured sets of variables to infer coordinate 
systems used for the neural control of movement has 
been a challenge. One of the dominant ideas originat¬ 
ing from the classical studies by Bernstein [5] has been 
that elements are united by the central nervous sys¬ 
tem into relatively stable groups to reduce the number 
of parameters manipulated at task-related neural lev¬ 
els. Such groups have been addressed as “synergies” [U 
138) or “modes” mm- Here, we are going to address 
such groups as “modes” to avoid confusion with an¬ 
other definition of synergies as neural structures that 
ensure stability of performance based on modes (re¬ 
viewed in usd- Some studies emphasized the relative 
invariance of the mode composition across tasks HUE 
I5B1 while other studies showed that the modes could be 
rearranged quickly with changes in the stability require¬ 
ments Earn . We believe that some of the disagreements 
may originate from using the same term for different 
coordinate systems. 

Within our scheme, measured variables produced by 
elements (e.g., digit forces, joint rotations, and mus¬ 
cle activations) are united into modes that are rela¬ 
tively stable across task variations. These modes may 
reflect preferred changes in the referent body configu¬ 
ration (cf. [10]) based on the person’s experience with 
everyday tasks. Mode composition is reflected in the 
structure of response to a noisy external input and can 
be reconstructed using matrix factorization techniques 
such as principal component analysis, factor analysis, 
independent component analysis, and non-negative ma¬ 


trix factorization (reviewed in [37]. Unlike many earlier 
studies, we do not assume that the number of modes 
(the dimensionality of the space where the control pro¬ 
cess takes place, aq in our notation) is smaller than the 
number of measured variables. It may be larger. For 
example, in our experiment, forces of four fingers were 
measured. The dimensionality of Xi may be higher cor¬ 
responding, for example, to the number of muscles or 
muscle compartments involved in the task. 

According to our main assumption, there is another 
coordinate system that allows ensuring stability of per¬ 
formance using local control organized about each axis. 
We suggest using a term “control coordinates” for this 
system. Unlike modes, control coordinates are sensitive 
to task changes, particularly to changes in conditions 
that affect stability of performance. When a person en¬ 
counters a novel task, however, he/she searches for an 
adequate set of control coordinates that would allow 
implementing local control rules. 

Our experiment showed that a quick reaction to un¬ 
expected perturbations acts along directions in the fin¬ 
ger force space that are close to the directions of finger 
force deviations produced by the perturbations. In con¬ 
trast, these reactions formed larger angles with vectors 
reflecting finger modes 0, eigenvectors in the space of 
finger forces that reflect finger force changes when a per¬ 
son tries to act with one finger at a time. This result 
corroborates the idea that a quick corrective actions 
are organized not along mode directions but along axes 
of another coordinate system, close to the ones along 
which the system shows the quickest deviation in re¬ 
sponse to the perturbation. 

4.2 Relations to the Uncontrolled Manifold and 
Referent Configuration Hypothesis 

Figure [4] offers a block diagram related to the control 
of the hand based on a few levels. At the upper level, 
the task is shared between the actions of the thumb 
and the opposing fingers represented as a single digit 
(virtual finger, Arbib et al. 1985) with the same me¬ 
chanical effect as the four fingers combined. Further, 
the virtual finger action is shared among the actual fin¬ 
gers (our experiments analyze four-finger coordination 
at that level). Even further, action of a finger is shared 
among a redundant set of muscles contributing to fin¬ 
ger’s action. And at the bottom level, each muscle is 
organized into a set of motor units by the tonic stretch 
reflex feedback that stabilizes the equilibrium point of 
the system “muscle + reflexes + external load”. Only 
the last level may be seen as based on relatively well- 
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known neurophysiological mechanisms with the thresh¬ 
old of the tonic stretch reflex (A) serving as a control 
(input) variable for the muscle [9]. 

In more intuitive terms, consider controlling the mo¬ 
tion of a donkey using a carrot. The carrot trajectory 
defines time evolution of referent coordinates for the 
head trajectory, to which the head is attracted. But the 
head cannot move without moving the legs. So, this big 
carrot is translated into a redundant set of smaller mini¬ 
carrots for individual legs; and then, into even more 
micro-carrots for the joints and muscles. Coordinates 
of such carrots for salient body parts represent referent 
body configuration. 

While the scheme in Figure [4] ensures some stability 
properties of the action, changes in the overall organi¬ 
zation of action (e.g., changes in the RC trajectories) 
may be needed if the task changes or there is a major 
change in the external force field. The general principles 
suggested in this paper offer a solution for the problem 
of stabilizing action in such conditions. 

A few recent studies have shown that, when a ma¬ 
jor change in the external conditions of task execution 
takes place, corrective actions are seen in both range 
(ORT) and self-motion (UCM) spaces with respect to 
task-specific performance variables Hi- Moreover, 
self-motion dominates, which, by definition, is unable 
to correct the action. These observations suggest that 
no single economy principle can form the foundation for 
such corrections. They allow interpretation within the 
set of principles suggested in this paper. Any pertur¬ 
bation is expected to induce large effects in less stable 
directions (those that span the UCM) as compared to 
more stable directions (ORT). According to the AMN- 
rule, corrective action is organized along the most nim¬ 
ble of the coordinates that allow local control. Since the 
projection of the “most nimble” coordinate to the UCM 
is expected to dominate, one can expect the corrective 
action to be primarily directed along the UCM as well. 


4.3 Reasonably Sloppy Control May be Good Enough 

Several recent publications presented arguments in fa¬ 
vor of the general idea that the CNS may not solve 
typical problems perfectly but rather use a set of sim¬ 
ple rules that lead to acceptable solutions for most 
problems [TM20l . Sometimes, the rules fail to solve spe¬ 
cific problems and then healthy people make mistakes, 
fall, mishandle objects, spill coffee, etc. We presented a 
particular instantiation of such a set of rules (based on 
the AMN-rule) and showed that these rules were able to 
stabilize action with high probability. The experimen¬ 
tal demonstration of relatively small angles between the 


vectors of perturbation-induced force changes and cor¬ 
rective changes in finger forces support the feasibility 
of the AMN-rule. 

Experimental studies of the effects of practice on 
stability of redundant systems have shown the existence 
of two stages (reviewed in [15]). First, when a person 
learns a new task, stability in relevant directions is de¬ 
veloped reflected in an increase in the relative amount 
of variance along directions that span the UCM for the 
salient performance variables. Then, when accuracy of 
performance reaches a certain ceiling, further practice 
leads to a drop of variance in those seemingly irrelevant 
directions. Why would a person stabilize directions that 
have no clear effect on overall task performance ? We 
addressed this issue in the section on perfectionism. 

Selection of a particular point (range) within the so¬ 
lution space has been discussed as resulting from opti¬ 
mizing the action with respect to some objective (cost) 
function [25]. Note that only one point on the solution 
hyper-surface is optimal with respect to any given cost 
function. Other points within the UCM violate the op¬ 
timality principle even though they lead to seemingly 
perfect performance. In a sense, large variance within 
the UCM combined with low variance orthogonal to the 
UCM implies that the person is accurate but sloppy. 

In the course of practice, when the person is as accu¬ 
rate as one can possibly be with respect to the explicit 
task, further practice may stabilize directions within 
the UCM to ensure that performance remains as close 
as possible to optimal with respect to a selected crite¬ 
rion. This is what we call “perfectionism”. Note that 
perfectionism is never absolute (“nobody is perfect”), 
and the system remains sloppy, but the degree of slop¬ 
piness can be reduced. 

We are grateful to Yen-Hsun Wu and Sasha Reschechtko 
for their help with the experiments. 
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